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We compare waveforms and orbital dynamics from the first long-term, fully non-linear, numerical 
simulations of a generic black-hole binary configuration with post-Newtonian predictions. The 
binary has mass ratio q ~ 0.8 with arbitrarily oriented spins of magnitude Si /m\ ~ 0.6 and 52 / m\ ~ 
0.4 and orbits 9 times prior to merger. The numerical simulation starts with an initial separation of 
r ?s 11M, with orbital parameters determined by initial 2.5PN and 3.5PN post-Newtonian evolutions 
of a quasi-circular binary with an initial separation of r = 50M. The resulting binaries have very 
little eccentricity according to the 2.5PN and 3.5PN systems, but show significant eccentricities of 
e ~ 0.01 — 0.02 and e ~ 0.002 — 0.005 in the respective numerical simulations, thus demonstrating 
that 3.5PN significantly reduces the eccentricity of the binary compared to 2.5PN. We perform 
three numerical evolutions from r m 11M with maximum resolutions of h = M/48, M/53.3, M/59.3, 
to verify numerical convergence. We observe a reasonably good agreement between the PN and 
numerical waveforms, with an overlap of nearly 99% for the first six cycles of the (I = 2, m — ±2) 
modes, 91% for the (I = 2,m = ±1) modes, and nearly 91% for the (I = 3, m = ±3) modes. The 
phase differences between numerical and post-Newtonian approximations appear to be independent 
of the (£, m) modes considered and relatively small for the first 3-4 orbits. An advantage of the 3.5 
PN model over the 2.5 PN one seems to be observed, which indicates that still higher PN order 
(perhaps even 4.0PN) may yield significantly better waveforms. In addition, we identify features in 
the waveforms likely related to precession and precession-induced eccentricity. 

PACS numbers: 04.25.Dm, 04.25. Nx, 04.30.Db, 04.70. Bw 
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I. INTRODUCTION 



The discoveries of quasars, AGN, and other black-hole 
driven astrophysical phenomena in the 1960's demon- 
strated that the most energetic astrophysical phenomena 
are powered by gravity in the strong-field regime. This, 
in turn, spurred a renewed interest in classical General 
Relativity. The second major milestone in the revival 
of the theory was the realization that when astrophys- 
ical black holes merge, they release incredible amounts 
of energy in the form of gravitational radiation, mak- 
ing them the brightest objects in the universe. During 
their last few orbits, merging black- hole binaries release 
energy with a peak luminosity of about 10 -3 c 5 /G, 10 23 
times the power output of the Sun. 

There are currently major experimental and theoreti- 
cal efforts underway to measure these gravitational wave 
signals. On the experimental side, these efforts required 
the construction of kilometers long interferometers, such 
as LIGO [l| and VIRGO @, sensitive enough to mea- 
sure arm length distance changes smaller than the radius 
of a proton. While on the theoretical side, these efforts 
required major advancements in signal extraction tech- 
niques and the theoretical modeling of the gravitational 
wave sources. Modeling the gravitational radiation from 
compact object sources has been particularly difficult, as 
they require solving the fully non-linear Einstein Equa- 
tions of General Relativity on powerful supercomputers. 
However, even with the rapid advancements in computer 



power, solving the two-body problem in General Relativ- 
ity proved to be remarkably difficult, requiring over thirty 
years of research for the field to mature. Then in 2005, 
two complementary and independent methods were dis- 
covered that allowed numerical relativists to finally solve 
the black-hole binary problem in full strong-field grav- 
ity EHi. 

The rapid progress and the number of new theoretical 
insights that followed these breakthroughs have trans- 
formed the field of numerical relativity (NR); turning it 
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and on our theoretical understan ding of black-binary 
spacetimes [H, [H, [H, [H, [H, Hi, H3, HE HI] . 

One of the breakthrough methods, the 'moving punc- 
ture' approach 0, H|, was adopted by a majority of the 
NR groups and has proven to be accurate for the neutron- 
star binary and mixed neutron-star — black-hole binary 
problems [6(1 HI| , as well as for black- hole configurations 
with more than two black holes [13, [63| . 

On the subject of black-hole binaries, the NR com- 
munity is in very good agreement concerning a variety 
of results. Black-hole binaries will radiate between 2% 
and 8% of their total mass and up to 40% of their angu- 
lar momenta, depending on the magnitude and direction 
of the spin components, during the last few orbits and 
merger [29l EH . |42| . [43|. In general, these binaries will 
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These recoils can be very large when the black holes in 
the binary have si gnif icant spin components in the orbital 
plane [30, 2Cil.T64j| (up to 4000 kms" 1 for astrophys- 
ical binaries [2Q| and even 10000 kms -1 for extremely 
close hyperbolic encounters |64J ) . which has astrophysi- 
cally important effects [H, [glJEfl, EH, H El]. The ob- 
servational consequences of these large recoil velocities is 
an active area of current research [27|, E3, El, Ei| E3, El] . 

Currently, one of the most important tasks of NR is to 
assist LIGO, VIRGO, and other interferometric observa- 
tories, in detecting gravitational radiation and extracting 
the physical parameters of the sources. Given the de- 
manding resources required to generate these black-hole- 
binary simulations, and the sheer volume of the seven- 
dimensional space of intrinsic parameters of black-hole 
binaries, we need to develop techniques to model ar- 
bitrary binary configuration based on numerical simu- 
lations in a carefully chosen sample of the parameters 
space, in combination with post-Newtonian and pertur- 
bative calculations. One of the most promising of these 
approaches involves determining the region of common 
validity of the numerical simulations and post-Newtonian 
expansions, with the goal of modeling the full wave- 
form using post-Newtonian waveforms for the initial in- 
spiral and numerical waveforms for the late-inspiral and 
merger. This method was pioneered with the use of the 
Lazarus waveforms [461 ] and has readily been pursued af- 
ter the breakthroughs in NR. 

Comparisons of numerical simulations with post- 
Newtonian ones have several benefits aside from the the- 
oretical verification of PN. From a practical point of 
view, one can try to parametrize deviations of the cur- 
rent 3.5PN expansions to fit the numerical results [r35l . 
H, E3, El, l69l. or directly propose a phenomenological 
description [701 ]. and thus make predictions in regions 
of the parameter space still not explored by numerical 
simulations. Another important application, from the 
theoretical point of view, is to have a calibration of the 
post-Newtonian error in the last stages of the binary 
merger. The first results of comparisons for equal mass, 
non-spinning binaries are encouraging [H. |49|. Pfll. l72l 1731] . 
Recently this analysis was applied to equal-mass, equal- 
spin binaries with the spins aligned with the orbital an- 
gular momentum (and thus non-precessing) [74l l75l [76| . 

In this paper we compare the numerical and post- 
Newtonian waveforms for the challenging problem of 
a generic black-hole binary, i.e. a binary with unequal 
masses and unequal, non-aligned, and precessing spins. 
The goal here is to evaluate accuracy of the current or- 
der of post-Newtonian expansions when including spins 
effects, as well as to develop new criteria for testing both 
numerical and post-Newtonian developments. 

The paper is organized as follows, in Sec. [IT] we review 
the numerical techniques used for the evolution of the 
black-hole binaries, in Sec. Mil we present results from the 



numerical evolution of two similar generic black-hole bi- 
naries, and in II VI we analyze and compare different wave- 
form modes as computed numerically and with the high- 
est available post-Newtonian approximation. Finally in 
Sec.[V]we present our conclusions. 



II. TECHNIQUES 

To compute the numerical initial data, we use the 
puncture approach [77j] along with the TwoPUNC- 
tures [78[ thorn. In this approach the 3-metric on the 
initial slice has the form 7^ = (ipBL + w) 4 <5 a b, where 
iJjbl is the Brill-Lindquist conformal factor, 5 a b is the 
Euclidean metric, and u is (at least) C 2 on the punc- 
tures. The Brill-Lindquist conformal factor is given by 
ipBL = 1 + Ei=i ^/(^l^ - ^il), where n is the total num- 
ber of 'punctures', is the mass parameter of puncture 
i (m^ is not the horizon mass associated with puncture i), 
and fi is the coordinate location of puncture i. We evolve 
these black-hole-binary data-sets using the LazEv [79] 
implementation of the moving puncture approach 
In our version of the moving puncture approach we re- 
place the BSSN [H, EH, EH conformal exponent (p, which 
has logarithmic singularities at the punctures, with the 
initially C 4 field \ — exp(— 4cf>). This new variable, along 
with the other BSSN variables, will remain finite pro- 
vided that one uses a suitable choice for the gauge. An 
alternative approach uses standard finite differencing of 
4> [|| . Recently Marronetti et al. [83] proposed the use of 
W = y/x as an evolution variable. For the runs presented 
here we use centered, eighth-order finite differencing in 
space [631 ] and an RK4 time integrator (note that we do 
not upwind the advection terms). 

We use the Carpet 84] mesh refinement driver to pro- 
vide a 'moving boxes' style mesh refinement. In this ap- 
proach refined grids of fixed size are arranged about the 
coordinate centers of both holes. The Carpet code then 
moves these fine grids about the computational domain 
by following the trajectories of the two black holes. 

We obtain accurate, convergent waveforms and horizon 
parameters by evolving this system in conjunction with 
a modified 1+log lapse and a modified Gamma-driver 
shift condition [4j, |85j, and an initial lapse a(t = 0) = 



2/(1 + V 



BL) 



The lapse and shift are evolved with 



(d t - (3 l di)a = -2aK, 
d t /3 a = B\ 
t B a = 3/4<9 f f Q 



j]B a . 



(la) 
(lb) 
(1c) 



These gauge conditions require careful treatment of x, 
the inverse of the three- metric conformal factor, near 
the puncture in order for the system to remain sta- 
ble 0, E^, H3]. In practice one sets a floor value for \ 
equal to one-tenth of its initial minimum value. This 
floor is only needed for the first ~ 5M of evolution. As 
shown in Ref . (86| , this choice of gauge leads to a strongly 
hyperbolic evolution system provided that the shift does 
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not become too large. In our tests, W showed better 
behavior at very early times (t < 10M) (i.e. did not re- 
quire any special treatment near the punctures) , but led 
to evolutions with lower effective resolution when com- 
pared to x- We chose r\ = 3 for the simulations presented 
here. 

We use AHFinderDirect (87j to locate apparent 
horizons. We measure the magnitude of the horizon 
spin using the Isolated Horizon algorithm detailed in [88| . 
This algorithm is based on finding an approximate rota- 
tional Killing vector (i.e. an approximate rotational sym- 
metry) on the horizon tp a . Given this approximate Killing 
vector (p a , the spin magnitude is 



S [v] - — 



Sir 



(<p a R b K ab )d 2 V, 



(2) 



AH 



where K ab is the extrinsic curvature of the 3D-slice, d 2 V 
is the natural volume element intrinsic to the horizon, 
and R a is the outward pointing unit vector normal to 
the horizon on the 3D-slice. We measure the direction of 
the spin by finding the coordinate line joining the poles 
of this Killing vector field using the technique introduced 
in [43] . Our algorithm for finding the poles of the Killing 
vector field has an accuracy of ~ 2° (see [43[ for details). 
Note that once we have the horizon spin, we can calculate 
the horizon mass via the Christodoulou formula 



H 



S 2 /(4r 



(3) 



where m; rr = y/ A/(16ir) and A is the surface area of the 
horizon. 

We also use an alternative quasi-local measurement of 
the spin and linear momentum of the individual black 
holes in the binary that is based on the coordinate ro- 
tation and translation vectors [26[ . In this approach the 
spin components of the horizon are given by 



1 



5 M = a~ / <t>l]R°K ab d'V, 

°K J AH 



(4) 



where <j>y* = 8ejS m kr m e^ k , and r m = x m — x™ is the 
coordinate displacement from the centroid of the hole, 
while the linear momentum is given by 

P W = T- I %fa&( K <* - K lab)d 2 V, (5) 

87r J AH 

where G« = 8\. 

We measure radiated energy, linear momentum, and 
angular momentum, in terms of ^4, using the formulae 
provided in Refs. [9(|. However, rather than using 
the full ipi, we decompose it into I and m modes and 
solve for the radiated linear momentum, d ropp ing terms 
with I > 5. The formulae in Refs. [89|, 120] are valid 
at r = oo. We obtain highly accurate values for these 
quantities by solving for them on spheres of finite radius 
(typically r/M = 50, 60, • • • , 100), fitting the results to 
a polynomial dependence in I — 1/r, and extrapolating 



to I = 0, Ht| . Each quantity Q has the radial depen- 
dence Q = Qo + IQi + 0(l 2 ), where Qo is the asymptotic 
value (the 0(1) error arises from the 0(1) error in rtpi). 
We perform both linear and quadratic fits of Q versus 
I, and take Qo from the quadratic fit as the final value 
with the differences between the linear and extrapolated 
Qo as a measure of the error in the extrapolations. We 
found that extrapolating the waveform itself to r = oo 
introduced phase errors due to uncertainties in the areal 
radius of the observers, as well as numerical noise. Thus 
when comparing PN to numerical waveforms, we use the 
waveform extracted at r = lOOAf. The extrapolations of 
the radiated quantities are far more robust. 

We convert the (£, m) modes of -04 into (£, m) modes 
of h = h+ — ih x by calculating the Fourier transform of 
each mode, dividing by — u> 2 (where u> is the Fourier fre- 
quency), setting the value of the resulting transform to 
zero inside some specified window — uj w < lu < u) w , as well 
as chopping off the transform at frequencies larger than 4 
times the quasi-normal frequency, and finally taking the 
inverse transform. By setting the transform to zero in 
this window, we remove the spurious constant and linear 
terms from h (we also remove spurious high-frequency 
noise from the waveform by truncating the transform at 
~ 4 times the quasi- normal frequency) . We confirm that 
the calculation is correct by taking two time-derivatives 
of the resulting h and measuring how much the result- 
ing function differs from the original tpi (See Fig. [5] in 
Sec. IIIip . We also use an alternative waveform compari- 
son, based on the modes of i/j^ rather than h, which does 
not require this transformation. 

We compute the eccentricities of the orbits using the 
techniques of [9l[ and introduce a second technique based 
on Newtonian trajectories. In [9l|, the eccentricity eo is 
defined as 



e D (t) = 



r(t) ~ r c (t) 
r c (t) 



(6) 



where r c is obtained by fitting r(t) to a low-order polyno- 
mial in t 1 / 2 . The actual eccentricity eo is the amplitude 
of the oscillations in the resulting en (t) . We also intro- 
duce a second measurement of eccentricity e r defined by 



.(*) = r(t) 2 r(t)/M. 



(7) 



Here too, the eccentricity e r is the amplitude of the 
oscillations in e r (t). This formula for the eccentricity, 
which is only accurate for e -C 1, arises from the Newto- 
nian formula for the orbital radius r(t) = y/ M /Vt 2 (\ + 
esin(fit)) +0(e 2 ). Note that in both cases, e(t) has si- 
nusoidal oscillations and secular decay. The ellipticity 
is the amplitude of the sinusoidal oscillations, while the 
secular decay affects the accuracy of the ellipticity cal- 
culation when its large. However, by differentiating r(t) 
twice with respect to t, the secular terms are suppressed. 
Eq. [7] can be modified with higher PN corrections to 
yield 



ecos(m) « [r(t) -r (t)]/{rn 2 



(8) 
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where 



ro(t) 



n 2 = ^r[l-(3- r))(M/r) + 0(M/rf] , (9) 
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(1751 + 588??) (M/r) 



, (io) 



r (i) = [252 - (1751 + 5&8rj)(M/r)} , (11) 

105 r 4 

and ro(i) is the zero-eccentricity inspiral trajectory. 



A. Initial Data 

To generate the initial data parameters, we used ran- 
dom values for the mass ratio and spins of the binary 
(the ranges for these parameters were chosen to make 
the evolution practical). We then calculated approxi- 
mate quasi-circular orbital parameters for a binary with 
these chosen parameters at an initial orbital separation 
of 50M and evolved using purely PN evolutions until 
the binary separation decreased to 11M. The goal was 
to produce very low eccentricity orbital parameters at 
r = 11M, as suggested in 91]. This technique is rather 
different from the technique in [53|, which used multiple 
numerical evolutions to determine quasi-circular orbital 
parameters. The initial binary configuration at r = 50M 
had q = mi I ma = 0.8, Si/m\ = (-0.2,-0.14,0.32), and 
Si/ml = (-0.09,0.48,0.35). As described in Sec. El we 
used both truncated 2.5PN equations of motion for spin- 
ning binaries, and equations of motions including 3.5PN 
corrections (without the HsiS 2 ,3PN term). Our PN evolu- 
tions use the ADM-TT gauge which is the one closest to 
the numerical quasi-isotropic coordinates (to help reduce 
possible gauge ambiguities) [H, [§4| . We denote the two 
resulting configuration by G2.5 and G3.5, respectively. 
We then used the PN momenta, spins, and particle lo- 
cations to construct the initial data for the numerical 
evolution. We fixed the puncture masses by requiring 
that the total ADM mass be IM and that the mass ratio 
of the two holes has the specified value. We renormal- 
ized the parameters to obtain an ADM mass of lAf in 
order to aid comparison of the two configurations and the 
analysis. 

The initial data parameter are summarized in Table [U 
We evolved these data using our eighth-order (in space) 
accurate code. We evolved the G2.5 configuration us- 
ing 12 levels of refinement, with a finest resolution of 
h = M/48, Af/53.33, and Af/59.33, and the outer bound- 
aries placed at 3072M. We used the standard 5th-order 
Kreiss-Oliger dissipation operator and six buffer zones at 
the refinement level boundaries. For the timestep, we 
chose a CFL factor of 0.5 for the inspiral phase, and then 
dropped the CFL by a factor of 0.95 during the merger 
phase. We reduced the CFL because otherwise the simu- 
lation proved to be unstable during the very fast plunge 
phase (due to a violation of the CFL stability condition 
for our evolution system). We evolved the G3.5 configu- 



TABLE I: Initial data parameters for the numerical evolu- 
tions. Parameters for configuration G2.5 were obtained from 
a truncated 2.5PN evolution of a binary starting with an or- 
bital separation of r = 50M, while parameters for config- 
uration G3.5 were obtained from an evolution with 3.5PN 
non-spinning corrections. The punctures have mass parame- 
ters mf , horizons masses (Christodoulou) mf , momenta ±p, 
spins Si, and both configurations have a total ADM mass 
A/adm- 



G2.5 



G3.5 



G2.5 



G3.5 



mJ/M 

mf jM 

xi/M 

Vi/M 

zi/M 

Sf/M 2 

Sf/M 2 

Sf/M 2 

p x /M 

10V/M 



0.40659 
0.44841 
3.32770 
-5.15410 
0.51835 
0.017896 
0.069204 
0.034786 
0.072919 
-5.4117 



0.40694 
0.44833 
-2.57272 
-5.57057 
-0.47758 
-0.036840 
-0.0050028 
0.069584 
0.080499 
-0.743105 



mJ/M 
mg/M 
x 2 /M 
V2/M 

Z 2 /M 

SUM 2 
SUM 2 
SUM 2 
p y /M 
Madm/M 



4.12328 
0.56054 
-2.66216 
4.12328 
-0.41468 
-0.066727 
-0.098217 
0.14722 
0.048074 
1.00000 



0.456072 
0.56106 
2.05867 
4.45696 
0.40645 

0.025826 
0.14951 
0.11050 

-0.036311 
1.00000 



ration with the same setup as the Af/53.3 G2.5 configu- 
ration, but chose an initial CFL factor of 0.475 (there was 
no evidence of any instability with this reduced factor). 



III. FULLY-NONLINEAR NUMERICAL 
WAVEFORMS AND TRAJECTORIES 

We calculated tpi using our original 4th-order accurate 
extraction code, and measured the convergence rate of 
the amplitude and phase of the waveform separately. In 
Fig. [1] we show the (£ = 2, m = 2) component of ip4 of 
the G2.5 configuration for the three resolutions. Note 
the excellent phase agreement until about t = 1400 M. 
The phase error increases exponentially during the last 
2 orbits. In Fig. [5] we show the convergence of the am- 
plitude of the (£ = 2, m = 2) mode. The amplitude 
shows between third- and fourth-order convergence as is 
apparent by rescaling the amplitude differences by 1.5098 
(fourth-order) and 1.35808 (third-order). As can be seen 
in Fig. [3l the phase error converges to eighth-order for 
t < 1200Af. Beyond t = 1200A/ (which is the begin- 
ning of the rapid plunge) the convergence falls to fourth- 
order, as is apparent from the rescaling of the phase dif- 
ferences by 2.30573 (eighth-order), 1.5098 (fourth-order), 
and 1.35808 (third-order). Note a convergence order of 
8 (up to t = 1200Af) implies that the error in the phase 
for the highest resolution run is less than 0.06 radians 
for t < 1200A/. In Fig. [4] we show the amplitude as a 
function of phase. The phase error in the waveform con- 
verges to higher order than the amplitude because it is 
sensitive to the phase error in the orbit, which, in turn, 
is a function of the convergence of the evolution code. 
The amplitude, however, appears to be sensitive to the 
extraction algorithm's numerical error. 

In Table HTl we show the radiated energy, angular mo- 



5 



0.002 
0.001 


-0.001 



-0.002 



V 4 (M/48.0) (A 




y 4 (M/53.3) (ft 




Qu /M/SQ ■'iA : T ; 










200 



400 



600 800 
t/M 



1000 



1200 



1400 



FIG. 1: The {£ = 2,m = 2) component of ^4 for the G2.5 
configuration for the three resolutions. Note the excellent 
phase agreement until about t = U00M. 
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FIG. 2: Convergence of the amplitude of the G2.5 (I = 2, m = 
2) component of ^4. The amplitude shows between 3rd- 
and 4th-order convergence (as demonstrated by multiplying 
the deviations in the amplitude by 1.358 and 1.5098, respec- 
tively) . 



mentum, and gravitational recoil versus resolution for the 
G2.5 configuration (sjj [90j. Here, extrapolation errors 
(to infinite radius) in the radiated energy and angular 
momenta dominate the finite-difference errors, while the 
extrapolation errors in the recoil appear to be similar to 
the finite difference errors. In particular V* ec has a no- 
ticeable finite-difference error. This can be understood in 
terms of the sensitivity of the out-of-plane recoil to the 
angle that the spin direction makes with the infall direc- 
tion at merger. Thus orbital phase errors in the plunge 
can lead to significant deviations in the out-of-plane re- 
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FIG. 3: Convergence of the G2.5 phase of the (t = 2, m = 2) 
component of ip4. The phase shows 8th-order convergence up 
to t = 1200M, decreasing to between 3rd- and 4th-order con- 
vergence during the plunge (as demonstrated by multiplying 
the phase deviations by 2.305, 1.5098, and 1.358 respectively). 
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FIG. 4: The amplitude of the G2.5 (£ = 2, m = 2) component 
of ip4, versus the phase. Note that the phase becomes more 
negative as t increases. 



coil [9fil96|. 

The radiated energy, angular momentum, and the re- 
coil velocity for the G3.5 configuration are given in Ta- 
ble IIIII The radiated energy and angular momenta are 
slightly larger for the G3.5 configuration than the G2.5 
configuration. Note that for both configurations, the ra- 
diated angular momenta in the x and y directions are too 
small to accurately measure. It should be pointed out 
that the quoted uncertainties in the radiated quantities 
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TABLE II: The radiated energy, angular momentum, and 
gravitational recoil versus resolution for the G2.5 configura- 
tion. The quoted uncertainties are due to extrapolation r — > 
oo. Note that this configuration has eccentricity er> ~ 0.02 
and e r ~ 0.01 



M/48 



M/53.3 



M/59.3 



S rad /M 

KecCkms" 1 ) 
KlcCkms- 1 ) 
KecCkms- 1 ) 



0.0512 ±0.0039 0.0513 ± 0.0036 0.0514 ± 0.0033 
0.018 ±0.021 0.017 ±0.020 0.014 ± 0.013 
-0.05 ±0.12 -0.05 ±0.12 -0.05 ±0.13 
0.4445 ± 0.0081 0.4478 ± 0.0103 0.4466 ± 0.0077 
-1.6 ±5.7 -6.9 ±6.0 -2.2 ±5.5 
78.36 ±6.51 75.75 ±2.95 71.47 ±0.24 
934 ± 31 1008 ± 24 947 ± 16 



TABLE III: The radiated energy, angular momentum, and 
gravitational recoil for the G3.5 configuration. The quoted 
uncertainties are due to extrapolation r — > oo. Note that this 
configuration has eccentricity en ~ 0.005 and e r ~ 0.002 



£ rad /M 
JIJM 2 

JrlJM 2 

KUkms" 1 ) 



0.0522 ± 0.0042 
-0.20 ±0.27 
0.051 ±0.057 
0.4551 ±0.0029 
26.3 ± 5.2 
103.0 ± 5.7 
1529.9 ±8.9 



for G3.5 are due to extrapolation to infinity. Additional 
uncertainties, due to truncation-errors are not included 
(although results from G2.5 indicate that the uncertain- 
ties in the radiated energy and angular momentum due 
to truncation errors are small compared to the errors due 
to extrapolation). 

As a final point, we show that our method for calculat- 
ing ft from tp4 using truncated Fourier transforms, yields 
a reasonable approximation to the original ip^ after differ- 
entiating twice. In Fig. [5] we show h and ^4 of the sub- 
leading (I — 2, to = 1) mode of the G2.5 configuration 
(See however the discussion concerning the amplitudes 
of ft in Sec. lIVBip . 



A. Eccentricity and Precession 

In Figs. [3 and[H we show the orbital trajectory for 
the G2.5 and G3.5 configurations. With the time direc- 
tion suppressed, we see excellent agreement between the 
trajectories at the three resolutions. This is similar to 
the excellent agreement in the amplitude versus phase 
of the (£ = 2, to = 2) mode. However, when includ- 
ing time, as can be seen in Fig. [51 there is a significant 
difference between the high and medium resolutions for 
t > 1200M. Also note in Fig. [9] the large eccentricity 
(apparent from the oscillations in r) for the G2.5 config- 
uration and that G3.5 has reduced, but still large, eccen- 
tricity. Thus, assuming that the PN series converges, we 
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FIG. 5: A comparison of h and ipA for the (I = 2, m = 1) 
mode for the G2.5 configuration. The plot demonstrates that 
the windowing procedure apparently does not contaminate 
the waveform to a significant degree. 



need to include still higher-order PN correction to obtain 
low-eccentricity initial data parameters. The reduced ec- 
centricity of G3.5 compared to G2.5, lends support to 
the hope that higher PN order will give low eccentric- 
ity data. Alternatively, to produce low-eccentricity data, 
one can try to use the iterative methods of [53|, which 
have been shown to work well for non-spinning binaries. 
Using the methods of [9l[ , we can calculate the eccentric- 
ity e£>(t), as shown in Fig. 1101 From the figure, we can 
see that the eccentricity of G2.5, which is er> ~ 0.02, is 
more than 3.5 times as large as the eccentricity of G3.5, 
which is ~ 0.005. Using the formula e r for the eccen- 
tricities (See Fig. EJ yields e r ~ 0.0088 for G2.5 and 
e r ~ 0.0022 for G3.5. However, as can be seen in the fig- 
ure, the eccentricity for G2.5 decays throughout the evo- 
lution, while the eccentricity of G3.5 (although smaller 
than G2.5) remains roughly constant for t > 600M . This 
is consistent with the results seen in Fig. [TBI which shows 
that the 3.5 PN prediction for the eccentricity does not 
decay with time for sufficiently close binaries and small 
eccentricities. In Ref. 9l|, they found that using PN pa- 
rameters from a PN-evolved inspiral (from r = 40M to 
r = 11M) reduced the eccentricity of the resulting binary 
from e = 0.01, for a quasi-circular binary at r — 11M, to 
e = 0.002. Here we see eccentricities after a PN-evolved 



inspiral to r = 11 M between 2.5 and 10 times as big. 

In Figs. [L21 and 1 131 we show r — x\ — x-i versus time for 
the G2.5 and G3.5 configurations after performing a con- 
stant rotation that maps the initial orbital motion onto 
the xy plane. Orbital plane precession drives the increase 
in amplitude of the z-component of r. The precession 
of the orbital plane is itself driven by the precession of 
the total spin of the binary. Thus we can measure the 
rate of orbital plane precession by looking at the com- 
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FIG. 6: The trajectory difference xi — X2 for the G2.5 con- 
figuration. Note the orbital plane precession and the very 
good agreement between trajectories at the different resolu- 
tions (the tracks from the different resolutions are not distin- 
guishable on this scale). 
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FIG. 7: An xy projection of the trajectory difference xi — X2 
for the G2.5 configuration. Note the very good agreement 
between trajectories at the different resolutions. The initial 
orbital plane is inclined with respect to the xy plane, making 
the orbit appear more eccentric. 



ponents of the black-hole spins as a function of time. 
In Fig. [T3] we show the components of the spin of the 
larger black hole as a function of time for the G3.5 con- 
figuration. Note that the precessional frequency is quite 
low, with the precession occurring on a timescale of order 
1000M; consistent with the time scale in the amplitude 
modulation of the rotated Z\ — Z2 trajectory component 



FIG. 8: An xy projection of the trajectory difference X\ — X2 
for the G3.5 configuration. The initial orbital plane is inclined 
with respect to the xy plane, making the orbit appear more 
eccentric. 
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FIG. 9: The coordinate distance r = |a?i — a?a| between punc- 
tures versus time for the G2.5 and G3.5 configurations. Note 
that the large eccentricity in the orbit (apparent in the os- 
cillation in r) is reduced by using the 3.5PN equations to 
generate the initial data. Unlike in Figs. [6] and [7] here the 
differences between resolution becomes apparent during the 
plunge. These differences drive the phase error. Also note 
that the G3.5 configuration merges more slowly. 
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FIG. 10: The eccentricity e D (t) of the G3.5 and G2.5 config- 
urations, as calculated using the techniques of [9lJ 
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FIG. 11: The eccentricity e r (t) of the G3.5 and G2.5 config- 
urations and the 3.5PN prediction for the G2.5 configuration 
(the 3.5PN prediction for G3.5 is a factor of 10 smaller than 
the NR prediction). The inset shows the 'eccentricity' at early 
times when gauge effects dominate the trajectories. Note that 
the eccentricity of G2.5 decays throughout the evolution while 
the smaller eccentricity for G3.5 remains roughly constant be- 
yond t ~ 630M. At later times the eccentricities of G2.5 and 
G3.5 begin to agree. 

in Figs. [T3il and [TBI Despite this long timescale, preces- 
sion can affect the waveform modes on shorter timescales 
via mode-mixing effects. That is, precession of the or- 
bital plane will cause our mode decomposition (which 
uses a fixed z-axis) to mix different modes (which oscil- 
late at different frequencies). This can lead to a beat- 
ing effect that produces amplitude oscillations visible in 
the waveform. In particular, when the orbital plane is 
aligned with the xy axis, the m modes have a frequency 
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FIG. 12: The coordinate displacement r — X\ — X2 between 
punctures versus time for the G2.5 configuration after per- 
forming a constant rotation that maps the initial orbital plane 
onto the xy plane. Precession is responsible for driving the 
amplitude of r z . 
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FIG. 13: The coordinate displacement r — xi — xi between 
punctures versus time for the G3.5 configuration after per- 
forming a constant rotation that maps the initial orbital plane 
onto the xy plane. Precession is responsible for driving the 
amplitude of r z . 



of ~ TOw or bit- Hence if the m = 2 and m = 1 or m = 3 
modes mix, the resulting system will have a beat fre- 
quency of ~ w or bit; the same frequency as that due to 
eccentricity. Thus, oscillations in the amplitude of the 
modes at the orbital frequency can arise both from pre- 
cession and ellipticity. We will come back to this point 
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FIG. 14: The components of the spin for the more massive 
black hole in configuration G3.5 as a function of time. The 
precession of the spin drive the orbital plane precession. Here 
the precession timescale is of order 1000M. 

in Sec. ITVBl 



IV. POST-NEWTONIAN EQUATIONS OF 
MOTION AND WAVEFORMS 

In order to calculate post-Newtonian (PN) waveforms, 
we need to calculate the orbital motion of the binaries. 
We use the ADM-TT gauge, which is the closest to our 
quasi-isotropic numerical initial data coordinates [93l.l94|. 
In this paper, we use two different approximate PN equa- 
tions of motion (EOM) based on [93, H \°$. To con- 
struct the EOM we use the Hamiltonian provided in [97| , 
with the additional terms provided injgl, [9!| , and the 
radiation- reaction force provided in 97]. We then use 
the standard techniques of the Hamiltonian formulation 
to construct EOM for the particle locations, momenta, 
and spins. In the first approximate EOM, we included 
the purely orbital Hamiltonian up to 2PN order, spin- 
orbit coupling up to 2.5PN order, and spin-spin coupling 
up to 2PN order (for the conservative part). That is to 
say, we use the Hamiltonian 

H R — iff), Newt + #0,1PN + -ffo,2PN 

+-ffsO,1.5PN + #SO,2.5PN + -ffsS,2PN ■ (12) 

Here we only include the leading order radiation reac- 
tion (dissipative) effect. We refer to the above EOM as 
the "truncated" 2.5PN EOM because there are terms up 
to 2.5PN order. For the second approximate EOM, we 
included the 3PN orbital Hamiltonian and 3PN spin(l)- 
spin(2) coupling in the ADM-TT gauge (99[, i.e., we use 



the Hamiltonian 

H F = H R + ifo,3PN + #SiS 2 ,3PN (13) 

(the £TsiS 2 ,3PN term was also computed in [lOOl . llOll Il02| 
in a different gauge). For the dissipative part, we added 
the 3.5PN (non-spinning) radiation reaction terms, as 
well as the leading spin-orbit and spin-spin coupling to 
the radiation reaction 97]. In the radi ation rea ction 
terms, we use the Taylor series of the flux |l03[|l04| ]. We 
refer to this second EOM as the 3.5PN EOM (in prac- 
tice the 3.5PN radiation reaction terms contribute to the 
orbital EOM at 6PN order). 

We then use the following procedure to construct hy- 
brid waveforms from the orbital motion. First we use 
the 1 PN accurate waveforms derived by Wagoner and 
Will 105] (WW waveforms) for a generic: orbit. By us- 
ing these waveforms, we can introduce effects due to ec- 
centricity and effects due to black-hole spins, including 
the precession o f the orbital plane. On the other hand, 
Blanchet et al. [l06[ ] recently obtained the 3PN wave- 
forms (B waveforms) for non-spinning circular orbits. We 
combine these two waveforms to produce a hybrid wave- 
form that includes the known higher-order corrections 
to the waveform. Note that, in the comparisons men- 
tioned below, the 'truncated 2.5PN' waveforms and the 
3.5 PN waveforms were constructed from the same WW 
and B expressions. Differences only arise because the 
'truncated' 2.5 PN waveforms are based on particle tra- 
jectories obtained from the 'truncated' 2.5 PN EOM. 

In order to combine the WW and B waveforms, we 
need to take into account differences in the definitions 
of polarization st ates and the angular coordinates (See 
Eqs. (73)-(75) of [ill for the definition of the WW po- 
larization states and Sec. 8 of [l06t ] for the definition of 
the B polarization states). The WW waveforms use the 
standard definition of GW polarization states, which are 
the same as those derived from the Weyl scalar, but the 
B waveforms use an alternate definition; leading to a dif- 
ference in sign for all the (£, m) mo des o f h. The angular 
coordinates in the B waveforms in [l06| ] are derived from 
circular orbits in the equatorial (xy) plane. To directly 
compare the NR and PN waveforms, we must add an in- 
clination to the B waveforms because in the generic case 
the orbital planes are inclined (with a time dependent in- 
clination angle) with respect to the a;?; plane. Hence we 
need to use the procedure developed in [73|, Il07l ] to trans- 
form the (£, m) modes of B waveforms into modes with 
respect to our rotated spin basis (We provide a simple 
derivation of these transformations in Appendix Ia"]) . The 
following is an outline of the procedure. Let L = f x p 
be the instantaneous orbital angular momentum, where 

L = i(sin6i cos$l, sinOi sin$i, cosOl), (14) 
r = r(sin © r cos <£> r , sin © r cos <I> r , cos © r ), (15) 
p = p(sin6 p cos<I>p,sm6 ) 3Cos<f > p,cos©p), (16) 

and L, r, p, ©l, $l, @ r , $ r , Q p , <j> p are functions of 
time. The first step is to rotate the orbital plane onto 
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the xy plane. Let R(a, 0, 7) be a general rotation defined 
by the Euler angles a, [3, and 7, where we first perform 
a rotation through angle a about the z axis, followed by 
a rotation through angle j3 about the y axis, and finally 
a rotation through angle 7 about the z axis (in practice, 
we never need to perform this final rotation). Thus a 
rotation R(— <I>l(£), — Qh{t), 0) transforms L and r into 
L' and f", where 

L' = L(0,0,1), (17) 
f = r(cos$ B (t),sin$ B (<),0). (18) 

The {i, m) modes of the B waveform, in a frame where 
the orbital plane is the xy plane, can be written in terms 
of cos$b(£)i sin $#(£), r, and w or bit- In order to calcu- 
late the (£,m) modes of h with respect to the numeri- 
cal coordinates (wher e the orbital plane is incl ined ), we 
use the results of [zl Il07| . As was shown in [107| . the 
spin- weighted s spherical harmonics in the numerical co- 
ordinates are related to those in the rotated coordinates 
(where the orbital plane is the xy plane) by 

F/ m (Q) = e^^e- WQ ^ m (-/3)e-^y/ m ,(fi'), 

m' 

where a, (3, and 7 are the rotation angles described above 
(note 7 = 0), and the phase factor e lsx arises from the 
transformation of spin- weighted function under a change 
of spin-basis. In [73| it was shown that K m s , m is indepen- 
dent of s (see Appendix |A1 for an alternative proof), and 
is thus given by [107| | 

K e m s , m (-P) = d e m , m (-0), (20) 

where d l m , m (0) is the Wigner d matrix given by, 

d l m , m {0) = y/(£ + m)\{t - m)\{t + m')\{£ - m')\ 

( 1 \h-\-m —m 

x y \~±1 

^ k\(l + m - k)\(£ - m' - k)\(m' -m + k)\ 



. 
am- 



2k-\-7n' — rn 





cos — 
2 



2i-2k-m' +m 



(21) 



where the sum over k is such that the factorials are non- 
negative. Since h = h'e~ 2lx , we have 



htm = j hY lm dQ. 

rn' 
rn' 

= J2^ m "* L dl, m (e L )h> em ,. (22) 

in ' 

The remaining complication arises from the fact that 
both the WW and B waveforms contain terms for a non- 
spinning circular orbit. To avoid adding the common 



terms twice, we subtract them from the B waveforms. 
First, using the 1PN WW formulae, we obtain the wave- 
forms from non-spinning circular orbits in the equatorial 
plane. We do this by applying the 3PN EOM for circular 
orbits to the WW waveform formulae. We then rewrite 
the waveforms in terms of the gauge invariant variable 
x, which is the normalized frequency. The B waveforms 
are given in terms of x, so we can identify those terms 
in the WW waveforms also present in B waveforms in 
a unique way. We then remove these terms from the B 
waveforms. For our generic case, we rotate the subtracted 
B waveforms modes and add them to the modes of the 
WW waveforms to obtain the hybrid waveform. Note 
that there are no significant gauge ambiguities arising 
from combining the WW and B waveforms in this way 
because at 1PN order the harmonic and ADM gauges 
are equivalent (and hence the WW waveforms are the 
same in the two gauges) and the B waveforms are given 
in terms of gauge invariant variables. 

Note that we calculate the spin contribution to the 
waveform through its effect on the orbital motion directly 
in the WW waveforms and indirectly in B waveforms 
through the inclination of the orbital plane. Other effects 
of spin and orbital plane precession on the waveforms are 
currently not known. 



A. Orbital motion and initial parameters 

Following the procedure detailed in [9l|, extended to 
spinning particles, we used purely post-Newtonian evolu- 
tions of a nearly quasi-circular binary with initial orbital 
separation r = 50M to obtain the positions, momenta, 
and spins for a non-eccentric binary with separation 
r ~ 11M. The idea behind this procedure is that one can 
specify quasi-circular parameters with very low eccentric- 
ity for binaries with large separations using the conserva- 
tive part of the Hamiltonian (i.e. solve for circular orbits). 
The subsequent PN evolution then provides the PN pa- 
rameters (including radial momentum) of a closer binary 
with similar (but lower) eccentricity. The initial quasi- 
circular binary configuration at r = 50M had PN param- 
eters q = mi/m 2 = 0.8, Si/ml = (-0.2,-0.14,0.32), 
and S 2 /m 2 2 = (-0.09,0.48,0.35). We refer to the binary 
configurations obtained using the truncated 2.5PN and 
3.5PN EOM as G2.5 and G3.5, respectively. It turns 
out that the order of the PN evolution is critical for pro- 
ducing low eccentricity binaries. The eccentricity of the 
G2.5 configuration, as measured by a subsequent 2.5PN 
evolution is quite small. However, both the numerical 
and 3.5PN simulations, show a that the eccentricity for 
G2.5 is actually relatively large. Similarly, the eccentric- 
ity of the G3.5 configuration, as determined from the full 
numerical simulation, while smaller than the G2.5 config- 
uration, is still relatively large. We used these r ~ HAf 
parameters in our numerical and subsequent PN evolu- 
tions. 

It is interesting to note that in the generic case, the 
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FIG. 15: e r versus radius for a binary with spins aligned with 
the angular momentum. Here the eccentricity decreases with 
r for all radii. 

eccentricity, according to 3.5PN does not decrease with 
time at smaller radii. To demonstrate this, we show the 
eccentricity, calculated using the formula e r (t) = r 2 r/M, 
where the magnitude of the oscillations in e r (t) is the ec- 
centricity. In Fig.[T5]we show the eccentricity versus time 
for a configuration with the same spin-magnitudes and 
mass ratio as our generic case, but with the spins aligned 
with the orbital angular momentum. As can be seen, the 
eccentricity decreases with radius. However, in Fig. [16] 
we show the eccentricity calculated for our configuration, 
and one slightly modified to give an even lower initial ec- 
centricity, versus time. Here we see that the eccentricity 
decreases to about e ~ 0.0005 and then remains constant. 
On the other hand, for the low-eccentricity data, the ec- 
centricity actually increases until reaching e ~ 0.0005. 
From the figures is apparent that precession affects in- 
duce an apparent ellipticity to the binary's motion that 
is not radiated away (at least to this order in the PN 
expansion) . 

For the G2.5 configuration we used a truncated 2.5PN 
evolution, which began at r = 50M, to obtain the 
PN parameters provided in Table IIVI The specific 
spins of the two holes are S\/m\ = 0.3945883931 and 
S 2 /ml = 0.6008327554, respectively. The 2.5PN ADM 
mass, Madm = mi + nT>2 + H R , for these parameters 
is Madm/M = 0.9925682736, where H R is given by 
Eq. ifU). 

When using these parameters in the numerical evo- 
lution, and subsequent PN evolutions starting from 
t/Madm = 11.08236108, we normalized the PN param- 
eters by the ADM Mass (i.e. we use the parameters 
r -» vj -Madm > P -» p/Madm, and S -> S/M£ DM ). This 
renormalization is helpful because we choose to normal- 
ize our numerical simulations such that the total ADM 
mass is 1. However, due to the spurious radiation on 
the initial slice, the numerical black-hole masses change 
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FIG. 16: e r versus radius for the G3.5 configuration and a 
very similar binary, with parameters chosen to reduce the 
(PN) initial eccentricity. Note that the eccentricity at r < 10 
is constant and roughly the same for both configurations. 



TABLE IV: PN orbital parameters for the G2.5 and G3.5 con- 
figuration at an orbital separation of r ~ 11M, as calculated 
directly from PN simulations starting at r = 50M. mi and 
m 2 denote the masses, x, y, and z denote the components of 
r — Xi — a?2, Pi (i = x,y,z) denotes the linear momentum, 
and Su and S 2 i denote the spin angular momenta. 





G2.5 


G3.5 


mi/M 


0.4455115640 


0.4455115640 


m 2 /M 


0.5568894551 


0.5568894551 


x/M 


5.9453450513 


-4.5976488271 


y/M 


-9.2084320770 


-9.9544694746 


z/M 


0.9260944396 


-0.8775891873 


Px/M 


0.0723766737 


0.0799120544 


Py/M 


0.0477169131 


-0.0360468994 


p z /M 


-0.0053715184 


-0.00073769138 


S lx /M 2 
S ly /M 2 
S lz /M 2 
S 2x /M 2 
S 2y /M 2 
S 2z /M 2 


0.0176308357 


-0.0365711851 


0.0681788517 


-0.0049664012 


0.0342713607 


0.0690768531 


-0.0657393278 


0.0256376428 


-0.0967624976 


0.1484228759 


0.1450366736 


0.1096979400 



with time, and eventually equilibrate to a mass ratio of 
q = 0.7993 (the uncertainty in the numerical masses of 
the two holes was 8m ~ 0.00003 at the highest resolu- 
tions). Thus in order to compare the PN and numerical 
waveforms, we need to account for this change in mass 
ratio. To do this, we modified our choices of mi and 
TO2 such that Mabm/M = 1 and q = mi/mi — 0.7993. 
However, because our two PN evolutions systems have 
different Hamiltonians, we needed to use slightly differ- 
ent values of mi/M^DM and /-Madm in each case. 
Note that the spin angular momentum is not affected by 
the spurious radiation to a significant level because the 
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spurious radiation is nearly axially symmetric about the 
two holes. For the truncated 2.5PN evolutions we used 



mi/Madm 
m 2 /M adm 



0.4486274928 . 
0.5612754821 . 



(23) 



i.e. from the equation Madm = 1 = (<7 + I)wi2 + 
H R (q, m.2), while for the 3.5PN evolutions we used 



mi /Madm 
to 2 /Madm 



0.4486635058 , 
0.5613205377, 



(24) 



i.e. Madm = 1 = (<Z + 1)^2 + H F (q,ni2). We verified 
that these changes in the masses have a negligible ef- 
fect on the eccentricity and waveforms according to the 
PN evolutions. We then used both the truncated 2.5PN 
and 3.5PN equations of motion to evolve this modified 
configuration from r rs 11M. We made one additional 
change in the truncated 2.5PN evolution of G2.5. In 
our original truncated 2.5PN evolution from r = 50M, 
we used a simpler form of the radiation reaction term 
based on PN expansion in the orbital parameters r and 
p. While in the subsequent evolution, we used a new ex- 
pression (consistent with the old expression to 2.5PN or- 
der in the Taylor expansion of the PN orbital parameters) 
based on an expansion in the orbital frequency [97f . How- 
ever, because we changed the EOM, the truncated 2.5PN 
evolution of the G2.5 configuration, which according to 
the original system had very-low eccentricity, now has a 
small residual eccentricity (see Fig. ITT]) . The radiation- 
reaction terms is directly related to the radial motion of 
the binary. Therefore, the radiation reaction force is very 
important to determine the quasi-circular configuration, 
and differences in the force have a strong effect on the 
motion. This is an indication that 2.5PN is not accurate 
enough to model the binary's motion in the r — 50M to 
r = 11M range. 

For the G3.5 configuration, we used a 3.5PN evolution 
(that did not include the HsiS 2 ,3PN term) from r = 50M 
to r = 11 M to obtain the orbital parameters provided 
in Table [BO The 3.5PN ADM mass of this system is 
Madm/M = 0.9927145092, and, once again, we renor- 
malized the PN parameters by the ADM mass. When 
evolving this system numerically, we used slightly altered 
values of the spin 



S lx /Ml DM 
Sly /M adm 

Si z /M AT)M 

S2x/M ADM 

S2y/M ADM 
*SWM ADM 



-0.0368395795, 

-0.0050028494, 

0.0695838052 , 

0.0258257964, 

0.1495121453, 

0.1105030087, 



(25) 



which introduced negligible changes in the waveforms 
and eccentricity. Here too, we find that the black holes 
absorb spurious radiation arising from the initial data 
that changes the mass ratio to 0.79937. To model 
this change in the 3.5PN evolution, we changed the mi 



and m 2 PN masses to mi /Madm = 0.4485829815 and 
m.2 /Madm = 0.5611706488. Here too, the changes to the 
masses do not affect the motion or eccentricity of the bi- 
nary according to the 3.5PN evolution. Thus, one s hould 
use an iterative procedure, like those in Refs. (49l. Il08l |. 
to reduce the eccentricity. 

According to the truncated 2.5PN evolution (with the 
new radiation reaction term based on the orbital fre- 
quency discussed above), the G2.5 configuration has a 
relatively small eccentricity, as is apparent in the small 
oscillations of the time dependence of the 2.5 PN orbital 
radius displayed in Fig. 1171 However, both a subsequent 
3.5PN evolution and the numerical evolution showed that 
these data were highly eccentric. In Fig. [T7]we see that 
both the 3.5PN and numerical simulations produce sim- 
ilar, large orbital radius oscillations (which are due to 
eccentricity). The G3.5 configurations, which has very- 
low eccentricity according to 3.5PN, as is apparent in 
the non-oscillatory behavior of the 3.5PN orbital radius 
seen in Fig. [THJ still shows relatively large oscillations 
in the orbital radius of the numerical simulation. Thus, 
using the 3.5PN equations of motion to generate low- 
eccentricity initial data reduces the eccentricity, but not 
nearly to the extent seen in non-spinning binaries [oit ] . 




FIG. 17: The evolution of the orbital radius for the G2.5 con- 
figuration from the numerical, 2.5PN, and 3.5PN simulations. 
The residual eccentricity in the 2.5PN evolution is due to our 
using a different 2.5PN radiation reaction term from that used 
in the original evolution beginning at r = 50Af . Note that 
both 3.5PN and the numerical simulation indicate that this 
configuration has relatively large eccentricity. 



B. Comparison of NR. and PN waveforms 

We produced both 3.5PN and 2.5 truncated PN wave- 
forms for the G2.5 configuration and 3.5PN waveforms 
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FIG. 18: The evolutions of the orbital radius for the G3.5 con- 
figuration from the numerical and 3.5PN simulations. Here 
the numerical simulations shows that the eccentricity was re- 
duced, but is still relatively large, while the 3.5PN evolution 
indicates that the binary is non-eccentric. 



for the G3.5 configuration. In Figs. fT9l and 120] we show 
the real part of the (t — 2, to = 2) mode of the strain 
h for G2.5 and G3.5 respectively. Note the reasonable 
agreement of the numerical and 3.5PN waveforms for 
700 M in both configurations. The differences between 
the PN and numerical waveforms are larger than the nu- 
merical waveform errors at this time. Also note that 
the 3.5PN waveform shows evidence of an early merger 
and has a higher frequency than the numerical wave- 
form, while 2.5PN waveform shows the opposite behav- 
ior. In Figs. [21] and [221 we show the real part of the 
(I = 2, to = 1) mode of h for G2.5 and G3.5 respectively. 
Again, the agreement is fairly good at earlier times and 
3.5PN is more accurate than 2.5PN. Also, note the inter- 
esting oscillatory behavior of the amplitude of the real 
part of the (£ = 2, m = 1) mode for both configura- 
tions. Here the amplitude (of the real part) oscillates 
at about the precessional frequency (see Fig. [14]). For 
the (I = 3, m = 3) mode, we obtained results similar to 
the (1 = 2, rn = 2) mode, as seen in Figs. [23] and |2"41 
However, for this mode, oscillations in the amplitude are 
more pronounced. 



1. Amplitudes 

We are concerned with exploring two different effects, 
eccentricity and precession. Long-term precessional ef- 
fects, which modulate the amplitude of the waveform 
over many cycles, are more readily apparent in h be- 
cause differentiating h twice (to obtain -04 ) suppresses 
low-frequency oscillations in comparison to higher fre- 



FIG. 19: The real part of the (£ = 2,m = 2) mode of h for the 
G2.5 configuration from the numerical, truncated 2.5PN, and 
3.5PN simulations. Note that the 3.5PN prediction is closer 
to the numerical waveform and that 3.5PN predicts an early 
merger while 2.5PN predicts a late merger (as is evident by 
the amplitude of the mode versus time). 




t/M 

FIG. 20: The real part of the (£ = 2,m = 2) mode of h for the 
G3.5 configuration from the numerical and 3.5PN simulations. 
Here too, 3.5PN predicts an early merger (as is evident by the 
amplitude of the mode versus time). 



quencies. As the binary inspirals, the frequency of the 
oscillations increases with the orbital frequency. Thus 
there is a large ramp-up in the amplitude of ip4 near 
merger. This can mask other effects as we observe be- 
low. On the other hand, the transformation from ip4 to 
h can induce both high-frequency and low-frequency dis- 
tortions in h (i.e. numerical errors due to the windowing 
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FIG. 21: The real part of the (£ = 2, m = 1) mode of h for the 
G2.5 configuration from the numerical, truncated 2.5PN, and 
3.5PN simulations. Note the precession induced modulation 
in the amplitude of the oscillations. 
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FIG. 23: The real part of the {I = 3, m = 3) mode of h for 
the G2.5 configuration from the numerical, truncated 2.5PN, 
and 3.5PN simulations. Note the relatively high-frequency 
oscillations in the amplitude (roughly corresponding to the 
orbital period). 
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FIG. 22: The real part of the (I = 2, m = 1) mode of h for the 
G3.5 configuration from the numerical and 3.5PN simulations. 
Note the precession induced modulation in the amplitude of 
the oscillations. 



FIG. 24: The real part of the {I = 3, m = 3) mode of h for 
the G3.5 configuration from the numerical and 3.5PN simu- 
lations. Note the relatively high-frequency oscillations in the 
amplitude (roughly corresponding to the orbital period). 



procedure in the Fourier transform). Thus it is advan- 
tageous to compare both ifj4 and h between the PN and 
numerical simulations. 

In order to analyze the behavior of the (£, m) modes 
of the waveform, we decompose the modes into ampli- 
tudes and phases. In Fig. [25] we show the amplitude 
of the (£ = 2, m = 2) mode of h for the G2.5 config- 
uration. Here the 2.5PN waveforms appear to capture 
the overall amplitude behavior to better accuracy, while 
the 3.5PN waveforms capture the oscillations in the am- 
plitude. These oscillations occur at roughly the orbital 



frequency and are due mainly to eccentricity and, to a 
lesser extent, precession. As discussed above, precession 
can induce an oscillation in the (£ = 2, m = 2) mode 
at the orbital frequency by mixing the (I = 2, m = 2) 
and {I = 2,m = ±1) modes (and since the m modes 
have frequency ~ |m|w or bit, where w rbit is the orbital 
frequency, the resulting modes will show a beating ef- 
fect at the orbital frequency). A similar plot for the 
G3.5 configuration, Fig. shows that 3.5PN predicts 
very small amplitude oscillations, which seem to confirm 
that the oscillations seen in G2.5 are mainly due to ec- 
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centricity. Note that in Fig. [26] the amplitude of the 
numerical (£ — 2,m = 2) mode oscillates at about the 
orbital frequency with a significantly larger amplitude 
than the 3.5PN prediction; indicating that these oscil- 
lations are due to eccentricity (which is consistent with 
the relatively large oscillations in the numerical orbital 
radius). Since the transformation from tp4 to h can in- 
duce artifacts into the waveforms, it is also important to 
compare the PN predictions for ipn with the numerical 
waveforms. In Figs. [57] and [55J we show the amplitude 
of the (£ = 2,m= 2) of V>4 for the G2.5 and G3.5 con- 
figurations respectively. Note that, for ip4, 3.5PN gives 
a clearly better fit to the G2.5 waveform than truncated 
2.5PN. Note also that the agreement between the 3.5PN 
and numerical ip4 appears to be significantly better than 
the agreement in h. Thus it appears that the windowing 
procedure has induced a very-low frequency mode into h 
that yielded a net change in the amplitude of the wave- 
form. 

The effects of precession become apparent in the sub- 
leading modes h (and to a lesser extent, in the sub- 
leading modes of ^4). However, numerical errors in the 
lower amplitude modes are also more pronounced. In 
Fig. l2!Handl3T)lwe show the amplitudes of the (I = 2, m = 
1) mode of h for the G2.5 and G3.5 configurations, re- 
spectively. Here both 2.5PN and 3.5PN capture the sec- 
ular behavior in the amplitude nicely. Unlike for the 
(£ = 2, m = 2) mode, here the PN amplitudes oscillate 
much more strongly than the numerical amplitudes for 
the G2.5 configuration, while 3.5PN seems to capture 
both the short (orbital frequency) timescale oscillations 
and the longer (precessional) frequency oscillation (un- 
til t ~ lOOOAf ) for the G3.5 configuration. The damp- 
ing of the numerical oscillations for the G2.5 configura- 
tion are likely a consequence of the windowing procedure 
(which acts as a high-frequency and low-frequency fil- 
ter), as a similar damping is not apparent in ^4 (See 
Figs . [30l and [32]) . Although the G3.5 configuration has 
very low eccentricity (according to 3.5PN), the effects of 
eccentricity can increase as the binary separation falls be- 
low 15 M (See Fig. [16]). This effect appears to be related 
to precession because the eccentricity of non-precessing 
binaries (See Fig. I15[) decreases uniformly with binary 
separation. In addition, mode-mixing effects may also 
be partially responsible for these oscillations in the am- 
plitude of the (£ = 2, rn = 1) mode at the orbital fre- 
quency. The secular oscillation in the amplitude of the 
{£ = 2, to = 1) mode matches the precessional frequency 
(See Figs. [TBI and [30]) . and is thus likely a direct conse- 
quence of precession (the amplitude of the (£ — 2, m = 1) 
mode contain s sign ificant contributions from the spins, 
see Eq. (3) in [l09| ). 

The {£ — 2, to = 1) mode of ip4, as seen in 
Figs. |3"T1 and I3"2l again shows that the 3.5PN waveforms 
are clearly more accurate than the truncated 2.5PN wave- 
forms. The agreement of the 3.5PN waveforms for the 
G2.5 configuration is remarkable. Note that the long- 
timescale oscillation seen in the {£ = 2, m = 1) mode of 



h, which is likely due to precession, is not apparent in 
ip4 of the G3.5 configuration. However, as this effect is 
smaller in G3.5 (as seen by comparing Figs. |3"T1 and |3"2"]) . 
it may be hidden in ip4 by the ramp-up in amplitude of 
if} 4 near merger. 

Finally, in Fig.l3"3"landl3~Hwe show the amplitudes of the 
(£ = 3, to = 3) mode of h for the G2.5 and G3.5 configura- 
tions, respectively. An interesting feature of these modes 
is that the late-time amplitude oscillations, which are 
roughly at the orbital frequency, increase with time, in- 
dicating that they are due to the precession-induced late- 
time eccentricity apparent in Fig. [TBI For the G2.5 config- 
uration, 3.5PN produces a remarkably good fit, captur- 
ing all oscillations in the amplitude until t ~ 1400M . On 
the other hand, 3.5PN does not capture the early-time 
oscillations in the G3.5 configuration. A possible expla- 
nation for this result is that, as seen in Figs. [T71 and [T51 
both 3.5PN and the numerical simulation show similar 
eccentricities for the G2.5 configuration, but 3.5PN shows 
much lower eccentricity for the G3.5 configuration. This 
eccentricity leads to the early-time oscillation in the am- 
plitude of the {£ — 3, to = 3) mode that are not captured 
by 3.5PN. However, as the binary evolves, the effects of 
precession-induced eccentricity in the PN EOM increase 
and eventually dominate. This causes the amplitude of 
the oscillations in the 3.5PN waveform to increase and 
eventually become larger than the numerical amplitude 
oscillations. In Figs. EH and [M] we show the amplitude 
of the {£ = 3, m = 3) mode of ip4 for the G2.5 and G3.5 
configurations. Here too 3.5PN gives a remarkably good 
estimation for the amplitude of the mode. Note that the 
orbital-frequency oscillations seen in Fig.[M]are not read- 
ily apparent in Fig. 1361 (even in the PN waveforms). This 
shows one advantage of analyzing h over ^4; eccentricity 
and precessional effects are more apparent in h. 

From the amplitudes of each mode, we see that pre- 
cession and eccentricity impart signatures on the modes 
of the waveform at the orbital frequency. However, the 
long-time oscillations in the amplitudes, here apparent 
only in the (I = 2, to = ±1) modes, seem to be due purely 
to precession, and occur at the precessional frequency. 



2. Phases 

In Figs. |3"7I and |3"51 we show the phase differences be- 
tween the 3.5PN and numerical waveforms for the (£ = 
2, to = 1), (I = 2, to = 2), and (£ = 3, to = 3) modes. 
In all cases we normalized the phase differences by divid- 
ing by £tt. Note that we renormalize by £tt, rather than 
to7t. If the orbital plane were to lie along the xy plane, 
or equivalently, we chose spherical coordinates such that 
the 9 = corresponds to direction of normal to the or- 
bital plane, then we would expect the {£, to) modes to 
have frequency lu » to Worbit, and an error in the orbital 
phase of 5$ or bit would lead to an error in the phase of 
the {£, to) modes of to c^orbit- However, in that case 
the {£ = 2, to = 1) mode would be very small. Conse- 
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FIG. 25: The amplitude of the (£ = 2, m = 2) mode of h for 
the G2.5 configuration from the numerical, truncated 2.5PN, 
and 3.5PN simulations. The oscillations in the amplitude are 
much more pronounced in the numerical and 3.5PN simula- 
tions, indicating that these oscillations are likely due to ec- 
centricity. 




FIG. 27: The amplitude of the (£ = 2,m = 2) mode of ip4 for 
the G2.5 configuration from the numerical, truncated 2.5PN, 
and 3.5PN simulations. Note the very good agreement be- 
tween the 3.5PN and numerical waveforms. 
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FIG. 28: The amplitude of the (£ = 2,m = 2) mode of lj)4, for 
the G3.5 configuration from the numerical and 3.5PN simu- 
lations. 



FIG. 26: The amplitude of the (t = 2, m = 2) mode of h for 
the G3.5 configuration from the numerical and 3.5PN simula- 
tions. The amplitude oscillations in the numerical waveform 
are much larger than those in the 3.5PN waveform, indicating 
that they are likely due to eccentricity 



quently, in our non-aligned spin basis, the (£ = 2, m = 1) 
mode is actually dominated by contributions from the 
(£ = 2, m = ±2) modes (of the aligned spin-basis). Thus, 
in our configurations, the {£ = 2,m = 1) mode has 
frequency 2w or bit and the error in the phase scales like 
2<5<£> or bit- Note that the renormalized phase differences 



are qualitatively independent of the mode. We therefore 
focus on the {£ = 2, m = 2) mode. In Fig. [35] we show 
the phase difference between the numerical, 2.5PN, and 
3.5PN (£ = 2, m = 2) mode of h for the G2.5 configu- 
ration. From the plot we see that the phase difference 
improves with the higher PN order and changes sign. It 
thus appears that still higher-order PN corrections may 
make the waveform phases agree. As seen in Figs. [5TJ [TT)1 
and [23l the truncated 2.5PN phase evolution is slower 
than that of the NR and 3.5PN, and thus its phase lags 
behind the other two. The 3.5 PN evolution merges too 
quickly (but is still closer to the numerical evolution) and 
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FIG. 29: The amplitude of the {£ = 2, rrc = 1) mode of h for 
the G2.5 configuration from the numerical, truncated 2.5PN, 
and 3.5PN simulations. The secular oscillation in the numer- 
ical amplitude occurs at roughly the precessional frequency. 
Here the shorter-timescale oscillations apparent in the PN 
waveforms are much smaller in the numerical waveform. 



thus its phase leads the numerical one. 



3. Matching 

In order to quantitatively compare the modes of the 
truncated 2.5PN and 3.5PN waveforms with the numer- 
ical waveforms we define the overlap, or matching crite- 
rion, for the real and imaginary parts of each mode as 
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where R em = Re(he m ), hm = Im(h tm ), and 



< f,g >= 



f(t)g(t)dt. 



(28) 



Hence, Mf m = Mf m = 1 indicates that the given PN and 
numerical mode agree. To compare PN and numerical 
waveforms, we need to determine the time translation 
8t between the numerical time and the corresponding 
point on the PN trajectory. That is to say, the time 
it takes for the signal to reach the extraction sphere 
(here r — lOOAf). We determine this time translation by 
finding the time translation near St = 100M that maxi- 
mizes the agreement of the early time waveforms in the 
(£ = 2, m = ±2), (I = 2, m = ±1), and (£ = 3,m = ±3) 
simultaneously. We find St ~ 112, in good agreement 



FIG. 30: The amplitude of the (I = 2, m = 1) mode of h for 
the G3.5 configuration from the numerical and 3.5PN simula- 
tions. The secular oscillation in the numerical amplitude oc- 
curs at roughly the precessional frequency. Here the shorter- 
timescale oscillations (corresponding roughly to the orbital 
period) are present in both waveforms with very similar am- 
plitudes. 
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FIG. 31: The amplitude of the (£ = 2,m = l) mode of ^4 for 
the G2.5 configuration from the numerical, truncated 2.5PN, 
and 3.5PN simulations. Note the very good agreement be- 
tween the 3.5PN and numerical waveforms. 



with the expectation for our observer at r = 100M. 
We also determine an alternate time translation, one full 
wavelength in the (£ = 2, m = 2) mode longer, that in- 
creases the matching of the (£ = 2,m = 2) mode over 
longer integration periods. On the other hand, this new 
time translation, St = 233, causes the (£ = 3) modes 
to be out of phase, leading to negative overlaps. Thus 
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FIG. 32: The amplitude of the (£ = 2,m = l) mode of ij} A for 
the G3.5 configuration from the numerical and 3.5PN simu- 
lations. 
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FIG. 33: The amplitude of the (I = 3, m = 3) mode of h for 
the G2.5 configuration from the numerical, truncated 2.5PN, 
and 3.5PN simulations. Note the very good agreement be- 
tween the 3.5PN and numerical waveforms. Also note that the 
short-timescale oscillations (orbital period) grow with time at 
later times, indicating that, at least at later times, they are 
due mainly to precession. 




FIG. 34: The amplitude of the (£ = 3, m = 3) mode of h for 
the G3.5 configuration from the numerical and 3.5PN simula- 
tions. Note that the short-timescale oscillation at later times 
grow with time, indicating that these later-time oscillations 
are due to precession. The early-time oscillations in the nu- 
merical waveform (at the same frequency) are likely due to 
eccentricity. 
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FIG. 35: The amplitude of the (£ = 3, m = 3) mode of ip4 for 
the G2.5 configuration from the numerical, truncated 2.5PN, 
and 3.5PN simulations. Note the very good agreement be- 
tween the 3.5PN and numerical waveforms. 



by looking at the (£ = 2) and (£ = 3) modes simul- 
taneously, we can reject this false match. The results 
of these matching studies are summarized in Tables Ivl 
and I VII As seen in the tables, the matching of the 3.5PN 
and numerical waveforms are significantly better than 
the matching of the 2.5PN and numerical waveforms for 
all modes. Similarly, all PN modes match the numeri- 
cal waveforms better over the shorter integration times. 
This is consistent with the qualitative agreements in the 



waveforms seen in Figs. |T9] [24j Note that the 3.5PN 
and numerical waveform matches for all modes are sig- 
nificantly better for the G3.5 configuration than the G2.5 
configuration for the longer t = 1000M integration time 
(the differences between the matches are most striking 
for (£ = 3, m = ±3) modes, where the matching is ~ 0.7 
for G3.5 and ~ 0.4 for G2.5). The only place where the 
matches for the G2.5 configuration are consistently bet- 



19 



0.0001 




180 330 480 630 780 930 1080 1230 1380 
t/M 

FIG. 36: The amplitude of the (£ — 3, m = 3) mode of tp4 for 
the G3.5 configuration from the numerical and 3.5PN simu- 
lations. 
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FIG. 37: The phase differences in h between the numerical 
and 3.5 PN simulations for the G2.5 configuration in the (I = 
2,m = 1), (£ = 2,m = 2), and (I = 3,m = 3) modes. We 
multiplied the phase differences in the modes by a factor of 
l/(£n). We divide by in, rather than mn, because the (£ = 
2,m = 1) mode is dominated by mode-mixing from the (£ = 
2,m — ±2) modes (see text for more details). The vertical 
lines shows the times when the (I = 2, m = 2) frequency is 
Mlo = 0.05 (t ~ 323M) and Mlo = 0.075 (t ~ 1075M). 



ter than the matches for the G3.5 configuration is the 
is the (£ = 2, m = ±1) modes for the shorter integra- 
tion times. Thus, it appears that the 3.5PN waveforms, 
in general, produce superior results for the more circular 
G3.5 configuration, which is likely due to the fact that 
the higher PN order B waveforms are accurate for quasi- 
circular, rather than eccentric, binaries. 
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FIG. 38: The phase differences in h between the numeri- 
cal and 3.5 PN simulations for the G3.5 configuration in the 
(£ = 2,m = 1), (,£ = 2,m = 2), and (£ = 3,m = 3) modes. 
We multiplied the phase differences in the modes by a fac- 
tor of l/(£n). Note that the normalized phase differences are 
qualitatively independent of the mode and arise from the or- 
bital phase error in the PN approximation. We divide by 
In, rather than mn, because the (£ = 2,m = 1) mode is 
dominated by mode-mixing from the (I = 2, m = ±2) modes 
(see text for more details). The vertical lines show the times 
when the {t = 2, m = 2) frequency is Mlo = 0.05 (t ~ 360M), 
Mlo = 0.075 (t ~ 1252M), and Mlo = 0.1 (t ~ 1493M). 



In Tables IVlTll Villi we show the matching of the modes 
of ip4 between 2.5PN, 3.5PN, and the numerical tp4. Here 
we find a better match when we use a slightly altered 
time offset. Note that matching is generally worse than 
that observed with h, especially for the longer integration 
times. This is consistent with the observation that the 
amplitude of "04 increases more rapidly in time than h 
(due to the effects of increasing frequency and the two 
time derivatives). Thus a matching of ^4 emphasizes the 
disagreement between the PN and numerical waveforms 
at later times. Interestingly, the matching in G3.5 is 
significantly better than G2.5 for the 1000M integration 
time, particularly in the (I = 3, m = 3) mode, where the 
matching between the 3.5PN and numerical ^4 is 65% 
for G3.5 and only 14% for G2.5. 



V. CONCLUSION 

We analyzed the first long-term generic waveform pro- 
duced by the merger of unequal mass, unequal spins, 
precessing black-hole binaries (a shorter simulation of 
this kind, which led to the discovery of the very large 
recoil configuration, was reported in [lEj]). We demon- 
strated eighth-order convergence of the waveform phase 
and fourth-order convergence of the amplitude (consis- 
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FIG. 39: The phase difference in the (£ = 2, m - 2) mode of 
h between the NR and PN waveforms for the G2.5 configura- 
tion. The vertical axis denotes the number of orbital rotations 
derived from GW cycle. Note that the normalized phase dif- 
ferences are qualitatively independent of the mode and arise 
from the orbital phase error in the PN approximation. 



tent with the order of accuracy of the extraction routine) 
in the numerical results. These waveforms clearly show 
the effects of eccentricity and precession on the ampli- 
tude in the sub-leading (I = 2,m= 1) and (£ = 3, m = 3) 
modes. In particular, analyzing the (£ = 2, m = 1) mode 
provides a way of detecting precessional effects in the 
observed waveforms. We have also found that there are 
two sources of eccentricity for a generic binary. Resid- 
ual eccentricity, due to a non-ideal choice of initial data 
parameters that tends to damp out as the binary separa- 
tion decreases, and precession-induced eccentricity that 
grows as the orbital separation falls below ~ 15M (this 
increase in eccentricity at later times is apparent in the 
(I = 3, m = 3) mode of h in both the PN and numerical 
waveforms) . 

We have compared these waveforms with the truncated 
2.5 post-Newtonian waveforms, as well as the waveforms 
with the non-spinning 3.0 PN conservative and 3.5 PN 
radiative corrections. We find a good initial agreement 
of waveforms for the first six cycles, with overlaps of over 
97% for the (£ = 2,m = ±2) modes, 90%-98% for the 
(£ = 2, m = ±1), and over 90% for the (£ = 3, m = ±3) 
modes. This provides a natural way to match numeri- 
cal waveforms to the post-Newtonian ones with a time 
translation (the same for all modes) motivated by the 
physical location of the observer (See Fig. for in- 
stance). The agreement degrades as we approach the 
more dynamical region of the late merger and plunge. 
The disagreement begins in a region where the numeri- 
cal waveform is still very accurate. Thus it appears that 
the disagreement is mainly due to errors introduced by 
truncating the PN series. Hence the overlap should be 



TABLE V: The overlap (matching) of the real and imagi- 
nary parts of the modes of h of the G2.5 configuration for the 
truncated 2.5PN and 3.5 PN waveforms and the numerical 
waveforms for various integration times and PN time trans- 
lation St. In all cases, we start the integration just after the 
numerical initial data (spurious radiation) pulse leaves the 
system. 



Integration Time 


600 


800 


1000 


Truncated 2.5PN (St = 112.2) 


Re (£ = 2, to = 2) 


0.789 


0.615 


0.365 


Re (£ = 2, to = 1) 


0.705 


0.501 


0.292 


Re (£ = 3, m = 3) 


0.596 


0.286 


-0.038 


3.5PN (St = 112.2) 


Re (£ = 2, to = 2) 


0.975 


0.922 


0.693 


Im (£ = 2,m= 2) 


0.976 


0.924 


0.723 


Re (I = 2, to = -2) 


0.975 


0.922 


0.693 


Im (£ = 2,m = -2) 


0.978 


0.926 


0.723 


Re (£ = 2, to = 1) 


0.982 


0.938 


0.687 


Im (i = 2,m= 1) 


0.977 


0.924 


0.699 


Re (£ = 2, to = -1) 


0.984 


0.939 


0.707 


Im (£ = 2,m= -1) 


0.980 


0.933 


0.711 


Re (£ = 3, to = 3) 


0.908 


0.794 


0.418 


Im (£ = 3, to = 3) 


0.916 


0.795 


0.435 


Re (I = 3, to = -3) 


0.909 


0.782 


0.403 


Im (I = 3,m= -3) 


0.912 


0.794 


0.426 


3.5PN (St = 233.3) 


Re (£ = 2, to = 2) 


0.928 


0.803 


0.746 


Re (I = 2, to = 1) 


0.918 


0.800 


0.774 


Re (£ = 3, to = 3) 


-0.850 


-0.602 


-0.492 



TABLE VI: The overlap of the real and imaginary parts of the 
modes of h of the G3.5 configuration for the 3.5 PN waveforms 
and the numerical waveforms. In all cases, we start the inte- 
gration just after the numerical initial data (junk radiation) 
pulse leaves the system. 



Integration Time 


600 


800 


1000 


3.5PN (St = 112.5) 


Re (£ = 2, to = 2) 


0.986 


0.964 


0.895 


Im (£ = 2, to = 2) 


0.987 


0.962 


0.900 


Re (I = 2, to = -2) 


0.986 


0.964 


0.895 


Im (£ = 2,m= -2) 


0.987 


0.962 


0.901 


Re (£ = 2, to = 1) 


0.904 


0.912 


0.843 


Im (£ = 2, to = 1) 


0.916 


0.901 


0.820 


Re (I = 2, to = -1) 


0.920 


0.908 


0.833 


Im (£ = 2, to = -1) 


0.917 


0.903 


0.816 


Re (£ = 3, to = 3) 


0.938 


0.891 


0.738 


Im (£ = 3, to = 3) 


0.919 


0.868 


0.721 


Re (£ = 3, to = -3) 


0.931 


0.880 


0.733 


Im (£ = 3, to = —3) 


0.906 


0.857 


0.721 



improved significantly by including 3.0 PN and higher- 
order conserva tive and radiative corrections, including 
spin terms IH, GM DM Gil EH ■ 

In fact, our results indicate that higher-order PN cor- 
rections to the orbital motion may further increase the 
accuracy of the PN waveforms. Although, the PN ex- 
pansion has not yet been shown to converge, we do find 
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TABLE VII: The overlap (matching) of the real and imagi- 
nary parts of the modes of ip4 of the G2.5 configuration for 
the truncated 2.5PN and 3.5 PN waveforms and the numeri- 
cal waveforms for various integration times with the PN time 
translation St = 106.5 for the truncated 2.5PN and St = 113.0 
for the 3.5PN. In all cases, we start the integration after 
t=180. The integration time means that the end of integra- 
tion is the same as that used in the overlap of h. 



Integration Time 




600 


800 


1000 


Truncated 2.5PN (St = 106.5) 


Re (£ = 2,ra = 2) 




0.900 


0.744 


0.435 


Im (I = 2,m= 2) 




0.898 


0.717 


0.469 


Re (£ = 2,m = 1) 




0.824 


0.654 


0.408 


Im (I = 2,m = 1) 




0.851 


0.675 


0.431 


Re (£ = 3, m = 3) 




0.767 


0.472 


0.00578 


Im (i = 3, m = 3) 




0.776 


0.477 


0.0102 


3.5PN (St = 113.0) 


Re (£ = 2,771 = 2) 




0.980 


0.909 


0.519 


Im (£ = 2,m = 2) 




0.984 


0.916 


0.563 


Re (£ = 2,m = 1) 




0.982 


0.936 


0.544 


Im (£ = 2,m= 1) 




0.976 


0.921 


0.594 


Re (£ = 3, 771 = 3) 




0.906 


0.759 


0.150 


Im (I = 3,m = 3) 




0.906 


0.754 


0.140 


TABLE VIII: The overlap of the real and imag 


;inary parts 


of the modes of tp4 of the G3.5 configuration for the 3.5 PN 


waveforms and the numerical waveforms with St 


= 113.5. In 


all cases, we start the inte; 


^ration after t= 


180. The integration 


time means that the end of intej 


^ration 


is the same as that 


used in the overlap of h. 










Integration Time 


600 




800 


1000 


3.5PN (St = 113.5) 


Re (£ = 2,771 = 2) 


0.981 




0.962 


0.860 


Im (£ = 2,m = 2) 


0.983 




0.958 


0.876 


Re (1 = 2,m = 1) 


0.882 




0.927 


0.850 


Im (£ = 2,m= 1) 


0.853 




0.893 


0.811 


Re (i = 3, 771 = 3) 


0.869 




0.841 


0.640 


Im (I = 3,m = 3) 


0.868 




0.834 


0.649 



remarkably better agreement in ip4 between the PN and 
numerical waveforms when moving from a 2.5PN EOM to 
a 3.5 EOM. This would appear to underscore the need for 
higher-order post-Newtonian calculations of both spin- 
orbit and spin-spin terms (especially in the EOM). Spin 
effects first appear at 1.5PN order, producing the spin- 
orbit hangup effect [2^, [4l|. Other spin effects, such as 
those due to precession, generate more subtle effects in 
the waveforms, and require higher-order PN corrections 
to accurately model (while subtle, these effects are also 
responsible for the very large kicks seen in spinning bi- 
naries with the spins oriented in the orbital plane). Our 
results seem to indicate that calculating these higher- 
order correction may prove to be invaluable for generat- 
ing waveform templates from generic black-binary con- 
figurations. 



Acknowledgments 

We thank E. Berti, A. Buonanno, A. Gopakumar, and 
R. Porto for careful review of the manuscript, and we 
thank B. Krishnan for suggesting the technique to cal- 
culate h from ip4. We thank the referees for their many 
helpful suggestions and insights, which have greatly im- 
proved the paper. We gratefully acknowledge NSF for fi- 
nancial support from grant PHY-0722315, PHY-0653303, 
PHY 0714388, and PHY 0722703; and NASA for finan- 
cial support from grant NASA 07-ATFP07-0158. Com- 
putational resources were provided by Lonestar cluster 
at TACC and by NewHorizons at RIT. 



APPENDIX A: TRANSFORMATION OF THE 
(£,m) MODES OF SPIN- WEIGHTED FIELDS 
UNDER ARBITRARY ROTATIONS. 

Here we consider the spin-weighted spherical harmon- 
ics in two different angular coordinate systems, (9, (j>) 
and (9', (f>'), related to each other by a simple rotation. 
For convenience, we will use f2 to denote the coordinates 
(9, <f>) and dft — sin9d9 d(j> to denote the volume element 
on the unit sphere. To construct spin-weighted functions, 
we need to define a null dyad q A on the unit sphere obey- 
ing q A q A 



and q q~A = 2 (indices are raised and 
lowered with the unit sph ere m etric). Here we will use 
q A = dg + %l sin 9 (see |l!2f for a review of the sub- 
ject). Any two choices for the dyad q A and q' A can differ 
by at most a phase factor, i.e. q' A — e lx q A . A spin- 
weight s field J transforms as J — * J' = e tsx J under 
this change in spin-basis. Of relevance here are the two 
dyads q A = do + ij sin 9 d<f, and q' A = d$> + ij sin 9' dp . 
The choice of q A fixes the 3 operator on spin-weighted 
fields. 

The spin- weighted spherical harmonics are constructed 
as follows [l07l |. 



where 



' (*-M)i / ( 

(£+\s\)\ 



-iy<5 s Y hn (n) 
A Y tm (n) 



if s > 
if s < 



(Al) 



3/ = d 8 f 
§./ = dgf - 



sin 

i 

sin( 



2<W - s/cotf 
;<W + sf cot 9 



(A2) 



for a function / of spin- weight s. In the £1' coordinates 
and the corresponding q lA spin basis. Eqs. (|A1[) and (|A2[) 
take on an identical form, but with the O coordinates re- 
placed with the f2' coordinates. Let J be a spin-weighted 
s field of arbitrary spin-weight that can be decomposed 
into spin-weighted spherical harmonics. That is, 



t=\s\ m=—l 



(A3) 
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We define a spin-zero potential j, such that 

oo I 

3 = E E UmYtm(ft), 

e=\s\ m=-e 

where 



(A4) 



■ _ T l (t-\s\)l 

and p = (— l) s if s > and p = 1 otherwise. Hence 
i s j if s > 



(A5) 



J 



gi s ij if s < o 



(A6) 



Under a change of spin basis, J — > J' — e tsx J but j 
j' = j. Thus 



y s j if s > o 

3'l s l j if s < 



(A7) 



and 



oo < 



J ' = E E 



^±^pf tm Yl s m m, (A8) 



where 



Thus 



(A9) 



_ . j (£+\s\)\ 

Jim- Jem ]l ( £ _\ a \yP> 



w _ , (j+M)j 



(A10) 



where 



(All) 



and hence we can determine how the modes of J mix 
under a rotation of the coordinates by looking at the 
modes of j. 

It was shown in jl07l | that the relationship between the 
spherical harmonic modes Yf, m (yt) and Y^ m (p,'), where 
the f2' coordinates are obtained from the Q, coordinates 
by a rotation described by the Eulcr angles a, 0, 7 in 
Sec. IIV1 is given by 

1 

Y tm (Q)= £ e^'^>&, m (-/9) 5^,(00, 

m' — — I 

(A12) 

and hence 



Jfn 



jY em (Q)dn 



jY em (n)dn' 
= j E e^ ma+m ^d e rn , m (-(3)Y; m ,(n')dn' 

J m'=-t 
I 

= E ^ m ' a+mi) dL m {-P)]', ml - (A13) 

Finally, using Eq. (|A10[) we get 

1 

Jlm= E e^' a+m ^dl lm (-P)J' lm , (A14) 

m'=-£ 

which is independent of the spin-weight of J. 
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